
load /Users/manojnair/projects/subamrine_cable_circulation/Shanghai-Amakusa.mat
% Hi Manoj,
% Here is one file from Shanghai-Amakusa.  It has original time and voltage data (timedata and voltdata),
% voltage time series with days with AP>20 changed to NaN (voltdata_corrected) and resampled weekly data (voltdata_rs).
% Thanks,
% Natalie

L = isnan(voltdata_corrected);

data_array = voltdata_corrected(~L);
time_array = timedata(~L);
%data_array = corrected_data_array;

b1 = min(time_array):14:max(time_array)+10; % define knots (in days)


if length(b1) > 1,
    
    sp=spline(b1,data_array(:)'/spline(b1,eye(length(b1)),time_array(:)')); % spline fit
    
end;

v=ppval(time_array,sp);


plot(time_array, data_array);
hold on;
set(gca,'FontSize',16);
plot(time_array, v,'r','LineWidth',2);
plot(week, voltdata_rs,'g','LineWidth',2)
legend('Unfiltered','Spline (2 weeks)','Weekly Resample');
xlabels('Volts');
xlabel('Years');
ylabel('Volts');
title('Shanghai-Amakusa');

%% remove solar tidal signals

t = time_array; % get a copy of the time array

%define solar and lunar periods
periods = [4.0000    4.8000    6.0000    8.0000   11.9672   12.0000   12.4210   12.6583   23.9345   24.0000   25.8910  (365.25/2) * 24 365.25 * 24];

%create a model
x=[cos(2*pi*t(:)/(periods(1)/24) ) sin(2*pi*t(:)/(periods(1)/24) )];

for j=2:length(periods)
    x = [x cos(2*pi*t(:)/(periods(j)/24) ) sin(2*pi*t(:)/(periods(j)/24) )];
end

%fit with the data

[s,stats] = robustfit(x, data_array);

% aplitude and phase
P = angle(complex(s(2:2:end),s(3:2:end)));
A = abs(complex(s(2:2:end),s(3:2:end)));

% create the synthetic time series



for i = 1:length(P),
    yy1(i,:) = A(i) * sin(2*pi*t*24/periods(i) + P(i)) ;
end;


% remove the modelled time series

corrected_data_array = data_array - s(1) + sum(yy);

%% HAW1

load '/Users/manojnair/projects/subamrine_cable_circulation/PTA.mat' HAW1S_voltdata timedata


L = isnan(HAW1S_voltdata);
data_array = HAW1S_voltdata(~L);


t = timedata(~L); % get a copy of the time array

%define solar and lunar periods
%periods = [4.0000    4.8000    6.0000    8.0000   11.9672   12.0000   12.4210   12.6583   23.9345   24.0000   25.8910  (365.25/2) * 24 365.25 * 24];
periods = [4.0000    4.8000    6.0000    8.0000   11.9672   12.0000   12.4210   12.6583   23.9345   24.0000   25.8910 ];
%periods = [4.0000    4.8000    6.0000    8.0000   11.9672   12.0000  24];


%create a model
m=[cos(2*pi*t(:)/(periods(1)/24) ) sin(2*pi*t(:)/(periods(1)/24) )];

for j=2:length(periods)
    m = [m cos(2*pi*t(:)/(periods(j)/24) ) sin(2*pi*t(:)/(periods(j)/24) )];
end

%fit with the data
[s,stats] = robustfit(m, data_array);

% amplitude and phase
P = angle(complex(s(3:2:end),s(2:2:end))); % 
A = abs(complex(s(3:2:end),s(2:2:end)));  % BUG found. I had used arctan (x/y) 


% remove unreliable tidal lines
p_re = stats.p(2:2:end);
p_im = stats.p(3:2:end);

L = p_re > 0.05 | p_im > 0.05;
A(L) = 0;

% create the synthetic time series

if exist('y'),
clear y
end;

for i = 1:length(P),
    y(i,:) = A(i) * sin(2*pi*t*24/periods(i) + P(i)) ;
end;


voltdata_corrected = data_array - sum(y);

% the Sq corrected voltage data is not sigificantly different from
% uncorrected data

%% Compare  HAW2 and HAW3 data

% HAW1 data
% load '/Users/manojnair/projects/subamrine_cable_circulation/PTA.mat' HAW1S_voltdata timedata
load '/Users/manojnair/projects/subamrine_cable_circulation/PTA.mat' HAW1N_voltdata timedata

L = isnan(HAW1N_voltdata);
data_array = HAW1N_voltdata(~L);


%data_array = voltdata_corrected;
time_array = timedata(~L);

b1 = min(time_array):30:max(time_array)+10; % define knots (in days)
if length(b1) > 1,
    
    sp=spline(b1,data_array(:)'/spline(b1,eye(length(b1)),time_array(:)')); % spline fit
    
end;
v=ppval(time_array,sp);


plot(time_array, v,'b*-');
hold on;



% HAW3 data

load  '/Users/manojnair/projects/subamrine_cable_circulation/MAK.mat' HAW3_voltdata timedata
% remove bad data
HAW3_voltdata(1:62290) = NaN;
L = isnan(HAW3_voltdata);
data_array3 = HAW3_voltdata(~L);

%data_array = voltdata_corrected;
time_array3 = timedata(~L);

b1 = min(time_array3):30:max(time_array3)+10; % define knots (in days)
if length(b1) > 1,
    
    sp=spline(b1,data_array3(:)'/spline(b1,eye(length(b1)),time_array3(:)')); % spline fit
    
end;
v3=ppval(time_array3,sp);


plot(time_array3, v3,'r*-')



% HAW3 data

load  '/Users/manojnair/projects/subamrine_cable_circulation/MAK.mat' HAW2_voltdata timedata
% Bad data for HAW2        

HAW2_voltdata(43006 : 43104) = NaN;
HAW2_voltdata(44498 : 44825) = NaN;
HAW2_voltdata(55837 : end) = NaN;
L = isnan(HAW2_voltdata);
data_array2 = -HAW2_voltdata(~L); % original data

%data_array = voltdata_corrected;
time_array2 = timedata(~L);

b1 = min(time_array2):30:max(time_array2)+10; % define knots (in days)
if length(b1) > 1,
    
    sp=spline(b1,data_array2(:)'/spline(b1,eye(length(b1)),time_array2(:)')); % spline fit
    
end;
v3=ppval(time_array2,sp);


plot(time_array2, v3,'k*-')


